# Set ylim on divergence
sweepso <- c(200,201,415,417,631,633,847,849,1064,1065)
sweeps <- c(415,417,631,633,847,849,1064,1065)
ylimit <- 0.3
for(ii in 1:length(sweeps)){
  i <- sweeps[ii]
  pdf(paste("sim",i,".distogram.pdf",sep=""),width=6,height=6,pointsize=12)
  real.distogram <- read.csv(paste("../sim",i,"/sim",i,".real.distogram",sep=""),header=FALSE)
  medianslin <- NULL
  lowerlin <- NULL
  upperlin <- NULL
  for(distclass in 1:300){
    data <- real.distogram[real.distogram[,1]==distclass,2]
    medianslin[distclass] <- median(data) 
    a <- mean(data)
    s <- sd(data)
    n <- length(data)
    error <- qnorm(0.975)*s/sqrt(n)
    lowerlin[distclass] <- a-error
    upperlin[distclass] <- a+error
  }    
  for(biop in 1:10){
    plot(0,0,t="n",ylim=c(0,ylimit),xlim=c(1,300),xaxt="n",xlab="Spatial distance",ylab="Divergence")
    title(main=paste("Sim =",i," Biopsies =",biop))
    lines(c(1:300),medianslin,lty=1,col="black")
    lines(c(1:300),lowerlin,lty=2,col="black")
    lines(c(1:300),upperlin,lty=2,col="black")
    #abline(h=divergence[divergence[,1]==7300,2], col="red")
    
    all.data <- NULL
    bcolors <- rainbow(10)
    for(r in 1:10){ 
      data <- read.csv(paste("../sim",i,"/sim",i,".biopsy.sampling.distogram.",r,".",biop,".10.log",sep=""),header=FALSE)
      data[,1] <- round(data[,1])
      points(data[,1],data[,2],col=bcolors[r])
    }    
  }
  dev.off()
}
